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The Hunga Tonga—Hunga Ha’apai (HTHH) volcanic eruption in January 2022 generated catastrophic tsunami 
and contends for the largest natural explosion in more than a century. The main island, Tongatapu, suffered 
waves up to 17 m, and Tofua Island suffered waves up to 45 m, comfortably placing HTHH in the “megatsunami” 
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league. We present a tsunami simulation of the Tongan Archipelago calibrated by field observations, drone, and 
satellite data. Our simulation emphasizes how the complex shallow bathymetry of the area acted as a low-ve- 
locity wave trap, capturing tsunami for more than 1 hour. Despite its size and long duration, few lives were lost. 
Simulation suggests that HTHH’s location relative to urban centers saved Tonga from a worse outcome. Whereas 
2022 seems to have been a lucky escape, other oceanic volcanoes have the capacity to spawn future tsunami at 
HTHH scale. Our simulation amplifies the state of understanding of volcanic explosion tsunami and provides a 


framework for assessment of future hazards. 


INTRODUCTION 

The Kingdom of Tonga comprises an archipelago of 169 islands sit- 
uated at the northern end of the Tonga-Kermadec Arc (Fig. 1). As 
the Earth's fastest-converging, most seismically active subduction 
boundary, this arc has the highest known density of submarine vol- 
canoes (1, 2). Sedimentary evidence on the east (and possibly west) 
coast of the main island of Tongatapu document numerous tsunami 
events in prehistory and perhaps as recently as the mid-15th century 
CE (3, 4), suggesting that large tsunami are not out of the ordinary. 
During an early 2015 eruption, the Hunga Tonga—Hunga Haʻapai 
(HTHH) submarine volcano (part of the 5-km-diameter submarine 
Hunga caldera) formed one of Earth's youngest volcanic islands (5, 
6). On 15 January 2022, the Hunga caldera explosively erupted at 
magnitudes far beyond the series of precursory Vulcanian or sub- 
Plinian eruptions that began on 20 December 2021. The 15 January 
eruption triggered shockwaves and pressure-coupled tsunami 
noticed as far away as the Caribbean. This event contends for the 
fiercest phreatoplinian volcanic eruption in more than a century, 
rivaling Krakatau in 1883, and displaying a VEI (volcanic explosiv- 
ity index) = 6.3 (7, 8). 

We build forward from several recent studies that audit the at- 
mospheric disturbances (9-16) and the pressure-coupled tsunami 
(17-19) induced by the explosion. Unlike these global-scale 
studies, we focus on the conventional tsunami generated by the ex- 
plosive displacement of seawater, within the Tonga Archipelago 
where the impacts were most severe. By using the high-resolution 
shallow water bathymetry model of Purkis et al. (20), we are able 
to simulate the passage of the tsunami waveforms through the ar- 
chipelago at much higher fidelity than previously achieved. 
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Combining pre- and post-satellite imagery (optical and radar), 
high-resolution digital elevation models (DEMs) with light detec- 
tion and ranging (LIDAR), and on-site field and drone surveys, 
we have assembled 118 tsunami runup measurements for the 15 
January 2022 HTHH explosive volcanic event spanning 10 islands 
of the Tonga Archipelago (Fig. 1). On the basis of Tsunami Squares 
simulations (21-24) and empirical pressure-distance-yield formulas 
derived from atomic bomb tests (25, 26), all of these runup mea- 
surements can, to a reasonable extent, be fit by using an isotropic 
explosive source model consisting of one or more discrete explo- 
sions not exceeding 15 megatons (Mt). Ancillary information 
from barometer and tide records, ear- and eyewitnesses, video 
records, and well-documented data on pressure-related glass 
damage further constrains the time sequence and yields. We 
propose at least five blasts, two of relatively small size around 4:00 
UTC (<0.1 Mt) and three increasingly larger ones at 4:06 (0.5 Mt), 
4:18 (4 Mt), and 4:56 (15 Mt). 

We consider our work to be of prescient value. First, it informs 
on the tsunami hazard posed by future eruptions of HTHH and 
other submarine volcanoes with shallow vents situated along 
oceanic arcs. Second, we contend that HTHH is an excellent 
natural laboratory to test hypotheses and models that can be de- 
ployed elsewhere to better understand similar eruptions and subse- 
quent tsunami as preserved in the geologic record. 


Diverse modes of volcanic tsunami 
Many tsunami-generating mechanisms exist during volcanic activ- 
ity (3, 27), including explosions, pyroclastic and debris flows, earth- 
quakes, flank and lava bench failures, caldera subsidence, and 
atmospheric waves. Unlike earthquake-induced tsunami, multiple 
tsunami can emanate from a single volcanic eruption, yielding 
waves of varied size, travelling at different velocities with different 
decay rates. Such diversity frustrates our current early warning 
systems and our modeling methods. 

According to the simplest theory, a pressure front of the form 
P(x-Vt) travelling at velocity V in a uniform ocean of depth H 
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Fig. 1. The study area spans the southern portion of the Tonga Archipelago in the South Pacific. (A) Location of model domain. (B) The Tongatapu and Haʻapai 
island groups situate on a series of carbonate platforms mapped by the 100-m isobath (cyan polygons). The HTHH island and Hunga submarine caldera lies to the west of 
these platforms in deep water (yellow arrow). (C) Tongatapu is the main island of Tonga and the site of its capital, Nuku‘alofa (city limits, black polygon). Here, tide gauges 
are stationed on the Queen Salote and Vuna Wharfs (red dots). The white star marks the position of Tsunami Rock, a 1600-ton erratic boulder in the vicinity of Fahefa 


village presumed deposited by prehistoric tsunami. 


couples to a water wave n(x,t) like 


P(x—Vt) H 
Py  (V*—gH) 0 


You can see that the water wave motion follows the pressure 
history in shape and velocity but is scaled by a velocity differential. 
These are called pressure-coupled water waves. Equation 1 holds for 
both the direct blast waves and acoustic gravity (AG) waves despite 
the forms of the pressure front differing in size, shape, decay rate, 
velocity, and duration. The colossal HTHH detonation spawned at 
least two types of pressure-coupled tsunami, plus a conventional 
tsunami (19, 28). 

The first type of pressure-coupled tsunami emitted by HTHH 
was driven by the shockwave of the direct blast (9, 12, 17). Blast 
waves in the air are compression driven and expand at velocities 
slightly faster than the speed of sound (340 m/s). The blast- 
coupled water waves mimic the expanding blast front both in 
speed and shape. The expanding blast pressure decays at least as 
fast as distance ~*’”, so the blast-coupled water wave attenuates 
quickly. These types of waves are included in our simulation 
because they aid in constraining explosive yield, but they are rela- 
tively tiny at distance. For example, our overpressure model (Eq. 
3) suggests that at 60-km distance, a 10 Mt blast (P = 0.3; 
psi = 2068 Pa) expanding over a 2000-m-deep ocean would 
couple to a water wave just 4 cm high. 

The second type of pressure-coupled tsunami is driven by the 
AG waves. These “tsunami in the sky” comprise low-pressure, 
long-duration, dispersive, gravity-driven atmospheric wave trains 
that continuously impart energy to the ocean (29-32). Similar to 


n(x,t) = 
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the direct blast case, the atmosphere-ocean coupling of AG waves 
is rarely efficient because of the velocity differential between atmo- 
spheric waves (270 to 300 m/s) and typical tsunami waves (100 to 
220 m/s). AG pressure waves decay more slowly with distance than 
blast pressures waves (~distance°° versus distance~*/”), so they 
could be detected across the globe. Also of note is that atmospheric 
pressure waves can run over land, creating waves in isolated water 
bodies where normal tsunami would not reach. AG-coupled water 
waves, 5 to 10 cm high, thousands of kilometers away were reported 
to have arrived hours earlier than later “conventional” eruption- 
generated tsunami waves (movie S1) (17, 19, 32-34). These pres- 
sure-coupled tsunami are a curiosity but of no real near-field, 
local concern. 

Lagging behind the HTHH pressure-coupled waves were the 
conventional tsunami generated by the explosive displacement of 
seawater. These travel on their own, guided by local y (€ H) 
speeds (35). Conventional tsunami were HTHH's primary destruc- 
tive arm and the primary focus of our inquiry. 


RESULTS 

Yields and timing of the HTHH tsunamigenic blast 
sequence 

By all accounts, the phreatoplinian HTHH eruption of 15 January 
2022 involved a series of eruptive and explosive events that spanned 
hours (7, 36). Eruptive events typically persist for several minutes, 
during which eruptive products are actively expelled from vents and 
great volcanic plumes are cast into the stratosphere (37, 38). Explo- 
sive events, by contrast, constitute a large explosion of steam 
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generated by phreatomagmatic interactions [ie., molten fuel 
coolant interactions (MFCI)]. MFCI explosions are sensed at dis- 
tance as “sonic booms." Because they distort the entire atmospheric 
column, both eruptive and explosive events may be responsible for 
the widely documented AG waves (9-14, 16, 17, 33). However, we 
contend that in the case of HTHH, the explosive or blast events 
created the bulk of tsunami (28, 34). 

While eruptive and explosive events go hand in hand in the 
HTHH sequence, the timings, spatial extent, and relative strengths 
of the two do not necessarily correlate well. In time, explosions last 
seconds, whereas eruptive events last minutes. In space, the loci of 
MFCI explosions need not necessarily be the same as the sites where 
eruptive products vented. Hence, the timing, spatial extent, and 
strength of the eruptive events deduced from satellite and ground 
observations (8, 15, 16, 39-42) may not jive with those same features 
of explosive events extracted from tsunami data. Our goal is to 
develop a plausible time history of explosive events that explain 
118 observed wave runup measurements and the two available 
tide gauge records in Nuku'alofa, as illuminated by tsunami simu- 
lation. Ancillary information comes from eye- and earwitnesses, 
window breaking reports, and the barometric pressure gauge 
record from the Queen Salote Wharf in the port (Fig. 2A). 

Earwitnesses in Nuku'alofa (Fig. 1C) report hearing two blasts 
just before 4:00 UTC. Additional eyewitnesses report the first 
plumes rising from HTHH at 4:06 UTC, reaching an altitude >20 
km by 4:11 UTC, and spreading umbrella-like over Nuku'alofa by 
4:25 UTC (43). We name these Blasts 1 and 2, but according to data 
recorded by the Nuku'alofa tide gauges (Fig. 2B), they produced no 
meaningful tsunami. Two increasingly loud booms were heard 
within the next 30 min; we name these Blasts 3 and 4. Damaging 
tsunami waves from this pair most certainly struck Nuku’alofa's 
west coast by 4:30 UTC and the waterfront of the city center by 
4:36. The two tide-gauge records show the first wave arrival from 
Blast 3 at 4:26 (Fig. 2B). Tsunami simulation fixes wave travel 
time from the HTHH volcano to the Vuna and Queen Salote 
Wharfs at 18 and 20 min, respectively. As far as tsunami are con- 
cerned, Blast 3 must have occurred at 4:06 UTC. Given the precision 
of the gauge data and the certainty in the position of the Hunga 
caldera relative to the gauges, there is no wiggle room here. Further- 
more, simulations show that a blast source at the volcano creates two 
or three separate arrivals at the tide gauges, separated by about 5 
min (Fig. 2B). These later arrivals represent the rebound of the 
blast cavity. We interpret the first two peaks in the tide gauges as 
a record of Blast 3. About 12 min later, three or more excursions 
of up to 1.25 m are seen in both gauges. We interpret those as 
being sourced from larger Blast 4 at 4:18 UTC. Using blast scaling 
laws based on surface, free-air chemical and atomic explosions (25, 
26) and our concepts of blast-generated tsunami, we estimate the 
explosive yield of Blasts 3 and 4 at 0.5 and 4 Mt, respectively. 

Although Blasts 3 and 4 as fixed above do an adequate job in 
matching the first 30 min of the tide gauge records, these two 
blast events alone are insufficiently powerful to reproduce the ob- 
served wave runup data. For a wide stretch along the north coast of 
Tongatapu waves ran up to 3 m above sea level (28), yet the tide 
gauges on the Queen Salote and Vuna Wharfs in Nuku'alofa 
show excursions less than 1.25 m relative to current tide. Tide 
gauges often read low compared to runup because they measure 
only wave height (potential energy), while runup reflects both 
wave height and wave velocity (kinetic energy). In addition, both 
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gauges are in compromised sheltered positions versus open 
coasts. Regardless, other widely spaced wave runup observations 
confirm that Blast 3 and 4 are lacking. A still bigger Blast 5 is needed. 

The timing of Blast 5 is uncertain, but it can be bracketed. First, a 
weather station and communication tower situated on a 13-m-high 
ridge on western Tongatapu was toppled by Blast 5 tsunami (44). 
The station logged every 10 min but transmitted only on the 
hour. It last broadcasted at 5:00 UTC. Allowing 10 min of 
tsunami propagation delay from HTHH to the tower, fixes Blast 5 
after 4:50 UTC. This timing is further confirmed by observations of 
people evacuating and later returning to the area (44). Second, and 
for reasons unknown, the tide gauge situated on Queen Salote 
Wharf stopped functioning at 5:24 UTC, a few minutes after the 
largest peak-to-peak value (>3 m) was recorded. The tide gauge 
on the Vuna Wharf, however, functioned throughout and recorded 
a third parcel of waves with peak-to-peak values >2.5 m starting at 
5:15 UTC and persisting until at least 6:00 UTC. Allowing for wave 
propagation delays, we propose that Blast 5 at 4:56 UTC spawned 
these waves. It is tantalizing to think that Blast 5 waves themselves 
knocked the Queen Salote Wharf tide gauge out of commission, but 
this cannot be determined unequivocally. 

How do we constrain the strength of Blast 5? There are at least 
two threads to pull here. One aspect is the reports of windows break- 
ing in Nuku'alofa. Breaking windows requires a quantifiable change 
in pressure, which, with overpressure versus distance formula (Eq. 
3), constrains the explosive yield. In the 1960's the Federal Aviation 
Authority collected a vast extent of sonic boom-related damage by 
flying supersonic jets over Oklahoma City six times a day for 26 
weeks (45, 46). Enough information was collected to construct 
probabilistic estimates of window glass breakage versus peak over- 
pressure (Fig. 3A). Reading from those results, 5% of windows 
would break under an overpressure of 0.15 psi, and 30% would 
break under an overpressure of a 0.35 psi. Anecdotal reports from 
Nuku'alofa suggest that 5 to 30% seems a reasonable span for glass 
breakage in the city. At 65-km distance, 0.15 and 0.35 psi corre- 
spond to explosive yields of 5 and 15 Mt, respectively (Fig. 3B). 
The barometric data recorded by the Tongan Met Office confirm 
that multiple pressure variations of up to 25 mbar (0.35 psi) oc- 
curred in the main city of Nuku'alofa (65 km from HTHH; 
Fig. 2A). These signals, starting at 4:45, 4:58, and 5:44 UTC, associ- 
ate with eruptive events as hot ash and gas were expelled from the 
caldera, resulting in long-period episodes of low pressure develop- 
ing beneath the rising plume (Fig. 2A). Those signals are too long of 
period to break glass in themselves, but they do show that HTHH 
had the capacity to issue pressure variations strong enough to break 
glass. It is worth stating here that the pressure measurements in 
Fig. 2 are collected and averaged over 60-s intervals. Explosions 
related to glass-breaking pressure changes last just seconds and 
therefore become smoothed over or lost in the Nuku'alofa barom- 
eter data. 

The most direct thread to pull is the 118 wave runup height mea- 
surements that we have derived from fieldwork and analysis of time- 
separated remote sensing imagery and topography. After running a 
number of simulations, we fixed the 4:56 UTC Blast 5 at 15 Mt. This 
seemed to offer the most reasonable fit to all the anchored observa- 
tions and collected measurements. 


3 of 14 


EZOZ ‘PI [dy uo S10'a0uaI0s' mmm: sdyy woy popeoumoq 


SCIENCE ADVANCES | RESEARCH ARTICLE 


HTHH Event #2 Event #3 
. [^N a 
r C O 
14.72/1015 & ^) C  2j—_|_—_+__,_ 
aan gl 


| 


Pressure (psi/mbar) 


14,29/985 


B3 waves B4 waves B5 waves 
i es 
B 04:26 04:38 05:15 Vuna gauge (10 s) 


co 
u 


— Vuna gauge (60 s) 
— QSW gauge (60 s) 


05:24 
QSW gauge 
fails 


= 
© 


0.5 


Height above current tide level (m) 
o 


-0.5 

-1.0 

-1.5 SS. Oa EA 
0-0 e. li el 


fi l i i 
05:30 06:00 06:30 07:00 07:30 08:00 
Time UTC 


Fig. 2. Time-averaged surface station data recorded by the Tongan Met Office 
in the main port of Nukuʻalofa on 15 January 2022. (A) Barometric pressure 
records three long-period, low-pressure phases starting at 4:45, 4:58, and 5:44 
UTC. We interpret these to be caused by the development of a low-pressure 
zone beneath the plume rising above HTHH. Averaged over 60 s. Intervals, these 
barometric data adequately capture these pressure drops that last approximately 
20 min. Pressure anomalies generated by the tsunamigenic explosions are too 
short lived to be recorded. To time these explosions, we rely on ear- and eyewit- 
ness accounts and the signal of the arriving tsunami waves at the two Nukuʻalofa 
tide gauges. The tide data (B) come from two gauges, one on the Queen Salote 
Wharf (QSW) and the second 1.8 km away on the Vuna Wharf. Situated 65 km away, 
the waves spawned at HTHH take 20 min to transit to Nuku'alofa. Two small blasts 
(Blasts 1 and 2) are reportedly heard in Nuku'alofa at around 4:00 UTC, but these 
evidently do not generate tsunami. Earwitness accounts of blasts at 4:06 and 04:18 
UTC correspond to parcels of waves arriving at the two tide gauges starting at 4:26 
and 4:38, respectively. A third wave parcel is recorded arriving at 5:15 UTC, deliv- 
ering the largest peak-to-peak variation in sea level (>3 m) recorded in the series, 
preempting the failure of the Queen Salote Wharf gauge at 5:24. We contend that 
this parcel was generated from Blast 5 at 4:56 UTC. The available data suggest that 
the explosive yields from Blasts 3, 4, and 5 were 0.5, 4, and 15 Mt, respectively. 


Tsunami simulation 
We simulated the tsunami parented by three explosive events of 
HTHH: Blast 3 at 4:06 UTC, Blast 4 at 4:18, and Blast 5 at 4:56 
[Fig. 4, movie S2, and online material by Ward (47)]. Blast 3, the 
weakest of this trio, had a yield of 0.5 Mt. By 4:10 UTC (Fig. 4A), 
the waves spawned by it reached the 100-m isobath of the two car- 
bonate platforms hosting the Tongatapu and Ha'apai island groups. 
The 4:21 UTC time step (Fig. 4B) informs that Blast 3 waves 
reached the coastline of Tofua in the north and the west coast of 
Tongatapu in the south. At this early stage, simulated wave 
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Fig. 3. Blast pressures, broken windows, and explosive yield. (A) In 1964, the 
U.S. Federal Aviation Authority (FAA) flew supersonic jet planes over Oklahoma 
City to evaluate the effect of sonic booms on people and property. From this ex- 
ercise, the FAA developed probabilistic assessments of the likelihood of window 
glass breaking for a given overpressure (46). Witnesses report that the blast occur- 
ring shortly before 05:00 UTC broke windows in Nukuʻalofa, situated 65 km from 
the HTHH volcano. The FAA data instruct that overpressures of 0.35 and 0.15 psi 
would break between 30 and 5% of windows, a span that sensibly brackets the 
damage in Nuku'alofa. Data formulated for overpressures at varying offset distance 
(Eq. 3) reveal that blasts with yields of 15 and 5 Mt, respectively, would deliver 
these pressures at 65 km (B). 


runups are modest (<3 m). As Blast 3 waves reach Tofua and Ton- 
gatapu, the larger waves emanating from 4.0-Mt Blast 4 radiate away 
from HTHH, such that at 4:21 UTC, waves from both blasts are 
in play. 

By 4:33 UTC (Fig. 4C), Blast 4 waves broadly impact all the 
islands in the model domain. Wave runups as high as 22 m are sim- 
ulated to strike the southern coast of Tofua, while waves ranging 
from 3 to 13 m reach the western coast of Tongatapu. The main 
city of Nuku'alofa, situated in a sheltered shallow lagoon on the 
northern coast of the island, has yet to experience meaningful 
waves. On this point, our simulation is in harmony with the data 
recorded by the tide gauges. These report sea-level perturbations 
of <0.5 m at 4:33 and must owe themselves to Blast 3. The 
western coastline of ‘Eua has received wave runups of up to 8 m 
by 4:33 UTC, and all the western islands of the Ha'apai Group 
have been impacted by the Blast 4 waves, although the <3-m 
runups are small by comparison. The 4:33 UTC time step empha- 
sizes the “wrap around" behavior of the waves in the south of the 
model domain. Note how the waves, initially radiating southward 
from HTHH, interact with the shallow bathymetry of Tongatapu 
and the complex of seamounts that lie 30 km to its west and shoal 
within 65 m of sea level (white box in Fig. 4C). These shallows 
refract the wavefronts to an easterly direction, a process reinforced 
by the presence of ‘Eua Island, that ultimately succeeds in redirect- 
ing the southward travelling waves to a northerly path toward the 
eastern sides of the Tongatapu and Ha'apai platforms, albeit with 
waves greatly attenuated by this circuitous route. 

By 5:00 UTC (Fig. 4D), 5 min after Blast 5 (15 Mt. at 4:56 UTC), 
its wavefront is already 30 km from HTHH and will imminently 
visit the western side of the Ha'apai platform. 

By 5:14 UTC (Fig. 4E), the largest predicted wave runups (44.9 
m) occur on the western coast of Tofua. Such extreme runups on 
Tofua are corroborated by scour marks observed on the flanks of 
this steep (and fortunately uninhabited) volcanic island, as corrob- 
orated by time-separated remote sensing, drone, and field observa- 
tions (Figs. 5 to 7). By now, the many islands on the west side of the 
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Fig. 4. Tsunami simulations for the repetitive blasts of HTHH. Wave propagation shown at T = 4:10 (A), 4:21 (B), 4:33 (C), 5:00 (D), 5:14 (E), and 5:25 UTC (F) and color 
scale depict maximum wave runup at the coast. These six time slices encompass the creation and propagation of waves from the three 15 January 2022 tsunamigenic 
explosions of HTHH. Highest runups in the model domain (45 m) are on the coastline of the northern island of Tofua that suffered the full unattenuated brunt of the 
tsunami. Receiving waves that have not crossed the shallow Tongatapu and Ha'apai carbonate platforms, the southwest coasts of the islands of Tongatapu and ‘Eua also 
receive runups >15 m high. Most other island locations, being protected by fringing reefs and shallow lagoons, fare better. Note how Tongatapu and ʻEua refract waves 
from a southerly to a northerly path (e.g., broken white arrows on frames 4:21 and 4:33 UTC), thereby projecting tsunami onto the (sheltered) eastern flanks of the 


Tongatapu and Haʻapai platforms. Animation as movie S2. 


Haʻapai platform receive Blast 5 waves, the largest impacts that these 
islands have so far suffered, but all <4 m high. Similarly, the west 
coasts of Tongatapu and ‘Eua have by now received their largest 
wave runups, simulated to be 16 and 12 m for the two islands, re- 
spectively. The pair of tide gauges on the Nuku'alofa waterfront 
receive their largest peak-to-peak values at this time (>2.5 m; 
Fig. 8), while the simulation predicts similar wave runups of 1.5 
to 2.5 m along the Nuku'alofa waterfront (Fig. 4E). 

By 5:25 UTC (Fig. 4F), only the Blast 5 wavefront is active in the 
simulation. It has now transited the Ha'apai platform and swept the 
islands situated on its eastern flank (Ha'ano, Foa, Lifuka, Uoleva, 
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and Uiha; Fig. 1 for locations). Between 2 and 8 m high, these 
runups are only marginally higher than those attributed to Blast 
4. To the south, the Blast 5 wavefront has bent eastward around 
Tongatapu to impact ‘ʻEua. Its west coast now receives the highest 
simulated runup at 17 m, just 1 m higher than the maximum 
achieved by smaller Blast 4. Notably, the Blast 5 wavefront, that 
has taken the faster, but longer, wrap around southerly route, 
now merges with waves from the same blast that traveled directly, 
but slowly, across the Ha'apai platform. The merged Blast 5 waves 
adopt the same northerly track, as those from Blasts 3 and 4 that 
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Fig. 5. Post-tsunami disturbances along the southern flank of Tofua Island 
based on synthetic aperture radar, DEM, and optical imaging analysis. (A) Lo- 
cation of Tofua 90 km north of the Hunga submarine caldera (HTHH). (B) Nadir 
view of the southern flanks of Tofua from Canadian Space Agency's (CSA) Radarsat 
Constellation Mission (RCM) C-band HH polarization spotlight SAR (1-m resolution, 
inverted in its backscatter) at 46° incidence with semitransparent Maxar WorldView 
(WV-02) visible wavelength data (resampled to 1 m) and superimposed DEM con- 
tours every 100 m from sea level (in purple), as well as the 40 m (above-sea-level) 
contour in orange. (C) Oblique view ray-traced using the co-registered synthetic 
aperture radar (SAR), with the 1-m ground sample distance (GSD) DEM, with the 
ocean view from WV-02 acquired 2 September 2022, and the inverted RCM-SAR 
acquired by the CSA in Spotlight (beam FSL30) mode on Aug. 24th, 2022. Pink 
to magenta correspond to possible disturbance areas (PDAs) identified in the sat- 
ellite data. These putative tsunami-related disturbances require field validation for 
confirmation but are suggestive of runups up to 45 m in this steeply sloping 
region. GEDI and ICESat-2 LIDAR topography has been used to evaluate the steep- 
est slopes in this central, southern portion of the Tofua volcanic island. RCM-SAR 
data courtesy of RCM image at 2022 courtesy of the government of Canada. 


traveled 1 hour earlier are the last vestiges of the tsunami in the 
model domain. 

Blast 5 offers an opportunity to examine the anatomy of the 
waves generated by a 15-Mt submarine volcanic explosion. Just 1 
min after the detonation, the simulation shows the blast-coupled 
wave radiating symmetrically away from HTHH (Fig. 9A). This 
wave is already 20 km from HTHH by 4:57 UTC and travelling 
close to the speed of sound but only 28 cm high and decaying 
swiftly. Following behind is the conventional tsunami generated 
by the explosive displacement of seawater. Immediately after the ex- 
plosion, the transient blast cavity that becomes the tsunami is 6 km 
across, forming a wave 85 m high on the north side of HTHH and 65 
m high to the south (Fig. 9B). Although the blast is isotropic, the 
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topography of the Hunga caldera is not. For this reason, the size 
of the conventional tsunami varies with azimuth. 

Tracking the peak wave height of the Blast 5 tsunami as it reaches 
the shallows surrounding Tongatapu emphasizes the effectiveness 
of bathymetric attenuation (transect X-X'; Fig. 9B). Having travelled 
55 km from HTHH, the 60-m wave that initially radiated south from 
the blast is 10 m high when it reaches the northern flank of the Ton- 
gatapu platform. As the wave crosses the fringing reef that veneers 
this flank, the height of the shoaling wave height nearly doubles (18 
m). Then, because of breaking, the wave drops to only 3 m. Attenu- 
ation continues as the wave transits the 12 km of shallow water that 
extends out from the north coast of Tongatapu. The wave that even- 
tually beaches on the island is <2 m high. By contrast, the Blast 5 
wave that visits the south coast of Tongatapu, a portion of the 
island unprotected by shallow water, is 6x higher (12 m; Fig. 9B). 

Last, as a reality check for the overall simulation, consider table 
S1, which lists the yield, energy release, and tsunami energy associ- 
ated with Blasts 3, 4, and 5. Ata minimum, for the calculation to be 
energetically acceptable, the energy of the tsunami produced must 
be less than the total energy lost in the blasts. Explosions in general 
are not very efficient sources of water waves because of the large pro- 
portion of energy lost to heat production, ejection of debris, turbu- 
lence, and seismic wave generation. Simulations here transfer just 5 
to 7.5% of the energy toward tsunami generation. The calculations 
made for acceleration-capping steps detailed in Materials and 
Methods are energetically acceptable. 


Accuracy of the tsunami simulation in space and time 

We evaluated the performance of the simulation in space against 
118 sites [Fig. 7 (A and B) for locations] distributed along 505 km 
of island coastline in the model domain (Fig. 7C). The performance 
of the simulation through time was judged from predicted and re- 
corded water levels at the Queen Salote and Vuna Wharfs in the 
main city of Nuku'alofa (Fig. 8). At the 91 field-observed sites, sim- 
ulated and measured runup heights are in good accordance, with a 
few notable exceptions. 

For the Tongatapu Group and its surroundings, field measured 
runups measured at Fafa Island (site IDs 98 and 99; Fig. 7B) are up 
to 5 m higher than predicted by the simulation. In addition, for 
waves reaching the north coast of Tongatapu (sites 100 to 109), 
the simulated values are a bit low versus the field data. This under- 
prediction intimates that the simulation overattenuates tsunami 
waves crossing the shallow shelf to the north of Tongatapu. In par- 
ticular, the frictional coefficient cg in Eq. 7 may be too large. On the 
west coast of Tongatapu, the situation reverses. Here, the simulation 
predicts higher wave runups than measured. Site IDs 44 through 55, 
for instance, show overpredictions up to 5 m in the extreme case 
(sites 49, 51, 53, and 54). Possibly, the simulated waves are insuffi- 
ciently attenuated before impacting this coastline because of inade- 
quate resolution of our near-shore bathymetry model. Note, 
though, that the west coast of Tongatapu experienced the largest 
waves in both the model and in reality. This deviation is acceptable 
considering the good agreement at all the other sites around the 
main island and adjacent ‘Eua (sites 40 to 43). 

For the Ha'apai Group, the simulation tends to underpredict 
field-measured wave heights on the central and east side of this plat- 
form by about 2.5 m (for instance, sites 13 to 36). On the west side of 
the Ha'apai Group, however, the simulation is in tune (e.g., sites 8 to 
10). These facts again point toward an overattenuation in the 
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Fig. 6. Drone imaged post-tsunami disturbances along the south-southwest flank of Tofua Island. (A) Maxar Worldview-02 view of SSW Tofua, with arrow indicating 
the area of DJI drone imaging (14 October 2022). (B) Oblique ray-traced perspective view of the SSW flanks of Tofua developed from DJI drone data as a structure-from- 
motion DEM by NASA Goddard. Spatial resolution of this DEM is 8 cm, with ortho-image mosaic of drone images projected atop the topography. Of the 1500 m of 
coastline imaged, approximately 300 m show signs of recent disturbance (in the DEM), consistent with the anticipated impact of the 15 January 2022 megatsunami. 
Drone images (C to F) emphasize disturbances ranging in height from 30 to 70 m, including landslides (C, D, and D), fresh debris fans, and debris chutes with transported 


trees (E and G). DJI drone data were acquired by S.J.C. 


simulation for waves crossing expanses of shallow water, west to east 
in this case. 

We supplemented the 91 field-visited sites with an additional 27 
sites where runups were derived from time-separated satellite and 
drone imagery with stereogrammetric topography and LIDAR. As 
detailed in Materials and Methods, observations from orbit are less 
precise than observations in the field and should be given less em- 
phasis, yet satellite data provide calibration points in the remote lo- 
cations unreachable during the field campaign. Satellite-derived 
runups around the Tongatapu Group accord well with the simula- 
tion (sites 77 to 83), reemphasizing the competence of the model in 
the south. Runups from satellite for the central and east side of the 
Ha'apai Group are higher than simulated, once more suggesting ex- 
cessive shallow water attenuation (sites 29 to 31 and 34 to 39). For 
islands on the west side of the Ha'apai Group (sites 3 to 10), the 
simulation and the satellite measurements match. The highest 
runups, observed or simulated, are located on the steep coasts of 
Tofua Island. Here, the simulation predicts that the Blast 5 waves 
had runups that reached 25- to 30-m elevation on the Island's 
southern coast (site 1) and 44 m on the west coast (site 2). Compar- 
ison of satellite and drone data before and after the HTHH eruption 
suggests vegetation scars at this height (Figs. 5 and 6). Although 
Tongatapu is closer to HTHH than Tofua (65 versus 90 km), the 
south coast of Tofua is unique in the model domain in not being 
mantled by a fringing reef or situated atop a shallow carbonate plat- 
form. Incoming tsunami would not be attenuated before running 
up the steep sides of this still-active volcanic island. Remote and sur- 
rounded by sea cliffs, Tofua is a national park and uninhabited. The 
high runups on Tofua radically contrast those on the other islands 
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that average <5 m. Further reconnaissance of Tofua in the aftermath 
of the 2022 tsunami is warranted using field and additional LIDAR- 
based methods, given its uniqueness. 

Simulation performance in time between 4:15 and 6:00 UTC, 
was assessed against the data collected by the Queen Salote and 
Vuna Wharfs tide gauges, both situated in Nuku'alofa on the 
north side of the main island of Tongatapu. The main channel of 
the Vuna Wharf gauge delivers data averaged over 60 s. Additional 
data streams deliver maximum and minimum water levels at 10-s 
intervals. These provide estimates of peak-to-peak variations 
(Fig. 8A). The tide gauge on the Queen Salote Wharf (Fig. 8B) 
lacks the 10-s data stream, but it did deliver data averaged over 60 
s. until the instrument failed at 5:24 UTC. 

Note that the two records correlate well in time, with waves ar- 
riving at Vuna Wharf 2 min earlier than at Queen Salote Wharf. The 
former is approximately 2 km closer to HTHH than the latter. The 
simulation correctly captures this 2-min offset. The simulation 
reveals that the first two or three observed wave peaks originate 
from Blast 3 starting at 04:26 UTC followed by two or three Blast 
4 waves starting at 4:38 UTC. These gradually decline until 5:15 
UTC, when the largest peak-to-peak deviations are seen in both 
the simulation and gauge data. These mark the arrival of Blast 5 
waves. Only in this period (5:15 to 5:35 UTC) of maximum commo- 
tion induced by the Blast 5 waves do the gauge and the simulation 
diverge. The simulation predicted water levels up to 1.5 m higher 
than recorded at the Vuna Wharf. By 5:40 UTC, the predicted 
Blast 5 waves realign with the gauge record. The limited data 
from the Queen Salote Wharf gauge provides further confirmation. 
Again, the successive increase in peak-to-peak deviation of sea level 
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Fig. 7. Comparison of simulated and observed wave runups at 118 sites. ID numbers relate the sites mapped in (A) and (B) to the comparisons of observed versus 
simulated data graphed in (C). In the maps, sites where runups were observed in the field are depicted by circles, and those observed using time-separated remote 
sensing are triangles. These symbols are color-coded to depict the discrepancy between observed runups and those predicted by tsunami simulation. In (C), purple and 
brown bars represent runup observations made in the field and made using remote sensing, respectively. Runups from the tsunami simulation are red bars. In each case, 
the bars span two estimates of observed and simulated runup at each site (Materials and Methods for details). 
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Fig. 8. Comparison of simulated and observed water levels. Sea level data 
taken from Vuna Wharf (A) and Queen Salote Wharf (B) in the Nukuʻalofa 
Lagoon (Tongatapu). The tidal signal has been removed and the current tide 
level set to 0 m at 4:15 UTC. Simulated water levels (red lines) mimic the gauge 
data through the passage of the Blast 3 and 4 waves, but slightly deviate for the 
period 5:15 to 5:35 UTC following the arrival of the Blast 5 waves, before falling 
back into step by 5:40 UTC. 
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induced by Blast 3 and Blast 4 waves is captured, as is their gradual 
dissipation. You can see a hint of Blast 5 waves arrive at the Queen 
Salote Wharf at 5:17 UTC, shortly before the station went offline. 

The tide gauge records and the field-measured runups reveal a 
conflict. The former show peak offsets <1.25 m, while all the field- 
measured runups in the area (sites 100 to 109) fall in the range of 2 
to 3 m. As suggested above, tide gauges often read low compared to 
runup because they measure only wave height (potential energy), 
while runup reflects both wave height and wave velocity (kinetic 
energy). Still, the simulation cannot match both. It compromised 
by slightly unpredicting the field-measured sites and slightly over- 
predicting the tide gauges between 5:15 and 5:35 UTC. 


DISCUSSION 

Planetary-wide observations of long-period atmospheric distur- 
bances that emanated from the HTHH eruption (VEI = 6.3) have 
delivered a good understanding of the volcano’s eruptive sequence 
(9-15, 17, 18). By contrast, the blasts embedded in that sequence 
and the tsunami that they spawned are poorly constrained. Unlike 
whole atmosphere disturbances, direct blast pressures and extreme 
tsunami waves cannot be observed at distance. Public-domain, 
near-field data recorded from Tonga are scant. We rely on record- 
ings from one barometric station and two tide gauges on 
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Fig. 9. Blast-coupled waves and conventional tsunami issued from HTHH. (A) Catches the simulation at 04:57 UTC, 1 min after Blast 5. Travelling close to the speed of 
sound, the tiny blast-coupled water wave is now 20 km from HTHH. Following behind is the conventional tsunami generated by the explosive displacement of seawater. 
Note how the 15-Mt blast is sufficient to expose dry seafloor (orange). This "bottoming out” diminishes the efficiency of tsunami generation. The X to X' transect delivers 
the cross section through the simulation shown in (B). Here, the red line charts the simulated peak wave height during the time that Blast 5 was active (4:56 to 5:15 UTC). 
Note that the height of the conventional tsunami as it departs the HTHH crater varies from 85 m in the north to 60 m in the south. Note too how the tsunami interacts with 
the shallow Tongatapu platform. The wave crossing the fringing reef outboard Tongatapu Island is 18 m high but is efficiently attenuated to only 3 mas it enters the broad 
lagoon behind the reef. Waves that eventually beach on the north coast of Tongatapu are <2 m high. Unprotected by a fringing reef, the wave visited upon the south coast 


of Tongatapu is 12 m high. 


Tongatapu, situated 65 km from HTHH. These data come with lim- 
itations. The barometric data average over 60-s intervals, too long to 
capture short-duration (a few seconds) blast pressure variations. 
The tide gauge data also average over 60 s and record activity at shel- 
tered locations in Nuku'alofa Lagoon. Waves there probably under- 
represent free-field sea-level motions (48). Still, by casting a wide 
net, we can fix times and yields of the blasts that generated 
tsunami. This relies on several threads: (i) the booms heard by ear- 
witnesses in Nuku'alofa; (ii) reports there of breaking windows; (iii) 
relationships between applied pressure and window breakage estab- 
lished by the U.S. Federal Aviation Authority (FAA); (iv) the tide 
gauge records; (v) runup values from widely ranged field, drone, 
and satellite measurements; and (vi) the fact that the blasts were 
surface, open-air explosions, as quantified by Eq. 3. 

Taking the data as a basket, we argue a sequence of three blasts. 
The first at 4:06 UTC with a yield of 0.5 Mt, the second at 4:18 UTC 
with 4 Mt, and the third at 4:56 UTC with a yield of 15 Mt. We rec- 
ognize that many mechanisms during a volcanic eruption generate 
water waves (3, 27). Likely, to some extent, many of these were in 
play during the HTHH sequence. Still, we demonstrate that three 
isotropic explosive wave sources adequately explain the existing 
tsunami data without recourse to other more complex mechanisms. 
Furthermore, we stress again that blast events differ from eruptive 
events and their timings may differ. For instance Wright et al. (16) 
interpret the Nuku'alofa barometric data to place events at 4:36, 
5:10, 5:51, and 8:46 UTC (their extended data figure 8). We interpret 
these several 20-min-long, low-pressure signals to represent erup- 
tive phases of the volcano, not explosive phases. The latter are too 
short lived to be recorded in that data. Furthermore, the timings 
from the work of Wright et al. (16) cannot be reconciled with the 
arrival times of tsunami waves along the Nuku'alofa waterfront, as 
reported by eyewitnesses and recorded by tide gauges there. Obser- 
vations of AG modes, not tsunami, are central to the work of Wright 
et al. (16), so the different interpretations are understandable. 
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However, it is relevant to note that our estimate of 15 Mt for 4:56 
UTC Blast 5 is bracketed by the 9 to 37 Mt derived by Astafyeva et al. 
(40) from ionosphere perturbations for the combined yield from the 
explosive and eruptive episodes of HTHH, although lower than the 
25 Mt estimated by Pakoksung et al. (28) in their near-field 
tsunami model. 

The timing of largest tsunami waves witnessed during the volca- 
no's activity, the toppling of the weather station on western Tonga- 
tapu (44), and the tide gauge excursions (Fig. 2), all fix Blast 5 to 
4:56 UTC. Curious that the largest explosive event of the sequence 
falls nearly one-half hour after the largest eruptive events, which oc- 
curred between 4:15 and 4:30 UTC (9, 12, 13, 37). Curious too is that 
the largest explosive event of the sequence has a negligible expres- 
sion in the global teleseismic records (12, 13, 41, 42). Perhaps Blast 5 
was not the same phreatomagmatic explosion-driven detonation as 
Blasts 3 and 4? Maybe part of Blast 5's 15 Mt of tsunami-generating 
energy came from non-explosive sources? A sudden structural col- 
lapse of the eruption column followed by submarine pyroclastic 
density currents could be sufficiently tsunamigenic to explain the 
data, while being less likely to be recorded seismically. It is up to 
other researchers to explore the plausibility of this. 

Another explanation might be that near—ocean surface phreato- 
magmatic explosions such as Blast 5 (and Blasts 3 and 4 for that 
matter) simply do not couple well to the solid Earth and thus do 
not show strong seismic signals. Instead, these types of explosions 
may send much of their energy upward or transfer a large propor- 
tion into tsunami or hydro-acoustic waves. 

The third explanation as to why Blast 5 failed to register in the 
far-field records is that we may have overestimated its yield. Two 
possibilities exist here. First, multibeam sonar bathymetric 
mapping in 2015 and again in 2016 revealed the Hunga caldera 
(with HTHH on its NW rim) to have a mean water depth of 
~150 m across its entire floor (49, 50). Surveys following the 2022 
eruption reveal the caldera to languish >1000 m below sea level 
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(Shane Cronin, pers. comm., 16 January 2023). It is logical to 
assume that Blasts 3, 4, and 5 were issued over ever-increasing 
water depth as the caldera collapsed. Our simulation fixed water 
depth at the volcano of about 150 m and Blast 5 tsunami “bottomed 
out” (Fig. 9), thus limiting its size. If, at the time of Blast 5, deeper 
water existed there, such “bottoming out" may not have occurred 
and a smaller yield would be needed to make the same-size 
tsunami. The second possibility for overestimating the yield of 
Blast 5 is that it was not a surface, open-air explosion to which 
Eq. 3 applies. An explosion deeper in the water column, for instance, 
might be more efficient than the surface explosions of Blasts 3 and 
4. That is, the yield needed to produce a given overpressure and 
tsunami size might be less for Blast 5 than Eq. 3 suggests. 

An intriguing finding from the simulation is the effect of blast 
yield on wave height at various locations. Consider the animation 
of our simulation (Supplementary Material) and the time steps ex- 
tracted from it (Fig. 4). Tsunami that traveled an unhindered path 
impinged on islands with wave heights that correlate to blast yield. 
At 65-km distance, for instance, the west coast of Tongatapu re- 
ceives wave runups of 8 to 12 m from 4-Mt Blast 4 but 15 to 17 
m waves from 15-Mt Blast 5. Runups on the south coast of Tofua, 
90 km from HTHH, are 17 to 23 m high from Blast 4 and 30 to 45 m 
from Blast 5. These results suggest that tripling blast yield approx- 
imately doubles the size of the tsunami unattenuated by a fringing 
reef. For many islands atop the Tongatapu and Ha'apai platforms, 
however, waves extensively attenuate as they cross expansive shallow 
water (e.g., Fig. 9B). After this journey, blast yield and wave height 
apparently decouple. Ha'ano, Foa, Lifuka, Uolva, and Uiha (Fig. 1B) 
mantle the southeastern flank of the Ha'apai platform. By the time 
the waves have traversed 30 to 50 km of platform to Uiha, for in- 
stance, Blast 4 and 5 wave runups are nearly identical (runups 
average 2 to 4 m) despite Blast 5 being >3x larger. Because 
shallow water attenuation is proportional to wave height squared 
(Eq. 7), larger-in does not mean larger-out. 

The simulation also highlights the diversity of wave heights 
across short distances. On the main island of Tongatapu, Blast 5 
waves beach on the southern coastline up to 17 m. Just a few kilo- 
meters to the southeast, the simulation predicts heights of only 10 m 
(Fig. 4E), and 10 km to the east, the Nuku'alofa tide gauges record 
only 1.5-m waves. On Eua, average wave height on the west coast is 
15 m but only 5 m on the east coast, 7 km away. Similarly, satellite 
and drone imaging of Tofua (Figs. 5 and 6) show substantial vari- 
ability of disturbance along the mapped coastlines, suggesting that 
tsunami risk varies on short length scales. 

Also captured in the simulation is how each blast makes more 
than one wave. The first wave associates with the direct outward 
push of the blast. Later ones associate with the rebound and oscil- 
lation of the blast cavity and span several minutes (28). Remaining 
separate when crossing deep water, these waves interact once 
impeded by bathymetry. Rising to within 65 m of sea level, a 
complex of seamounts lies 30 km off the west coast of Tongatapu 
(labeled in Fig. 4C). Here, multiple Blast 5 waves coalesce as the 
leading front stalls over this topography. The seamounts then act 
like secondary point sources of waves. Susceptible to bathymetry, 
waves spawned by the blasts refract as they pass the aforementioned 
seamounts and bend around the south coast of Tongatapu. The to- 
pography conspires to redirect the southward travelling wavefronts 
onto a northerly path, ultimately leading the tsunami to wrap 
around the Tongatapu platform to the point where even the 
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sheltered eastern side of the Ha'apai platform is illuminated. 
Another way that waves reach the east side of the Ha'apai platform 
is via the more direct (but slower) route across the platform top. 
With sequential wavefronts emanating from the successive HTHH 
explosions but traversing the archipelago via different routes at dif- 
ferent speeds, Blast 4 and Blast 5 waves interact in the west of the 
model domain, although the former generated 40 min earlier. 
These results highlight how a single tsunami can remain “captured” 
by an archipelago, how waves from multiple blasts interact even 
when separated by hours, and how wrap-around behavior can 
result in sizable tsunami beaching in areas where they would not 
necessarily be expected. 

The 15 January 2022 eruption of HTHH (within the submarine 
Hunga caldera) holds lessons for both past and future tsunamigenic 
events along the Tonga-Kermadec arc. Tonga makes a compelling 
case study, as substantial evidence exists for tsunami inundating 
large areas of Tongatapu in precolonial history (3, 4). Considering 
sedimentological and archaeological evidence, Lavigne et al. (3) pin 
a major tsunami in the mid-15th century with runup heights up to 
30 m. The powerful Tu'i Tonga kingdom that then existed was se- 
verely affected, and an oral history remains more than 500 years 
later. We believe the 2022 HTHH tsunami to be of similar scale. 
Simulated and field-measured wave runups on the west coast of 
Tongatapu approach 20 m. On the western coast of Tofua, they 
are over 40 m. Our simulated peak height of the Blast 5 wave at 
HTHH was 85 m (Fig. 9), in good accordance with the 90 m pre- 
dicted by Heidarzadeh et al. (34). Fortunately, the urban center of 
Tonga, Nuku'alofa, situates in a sheltered lagoon. Had the popula- 
tion been dispersed along the exposed ocean-facing coastlines of the 
archipelago, a much higher death toll might be anticipated than the 
estimated six lives lost. 

The wave runups from the 2022 HTHH explosive event comfort- 
ably meet the criteria for a megatsunami and contend for the largest 
event anywhere in the past 100 years. For comparison, field records 
of the 1883 Krakatau (Indonesia) and 2011 Tohoku (Japan) tsuna- 
mis both indicate wave heights of 35 to 40 m (51-53), placing these 
events and HTHH ina similar league. The tsunami spawned by the 
2018 eruption of the Anak Krakatoa produced maximum wave 
runups of “only” 13 m (but killed 437 people) (54-56). Another 
megatsunami generated by the 1957 Andreanof Islands earthquake 
in Alaska (57) had wave runups in the range of 9 to 23 m (58), about 
half the maximum of HTHH. Likely the deadliest in recorded 
history, the 2004 Sumatra tsunami produced wave runups generally 
less than 10 m (59) and locally up to 20 m (60), killing 227,898 
people in 14 countries (61). 

With once-in-a-century wave runups, some exceeding 40 m, the 
death toll from HTHH would be expected to be far greater. The 
main factors that led to this, we suggest, are the quirk of the location, 
a worldwide pandemic, and increased evaluation drills and aware- 
ness efforts carried out in the years prior. The western coastlines of 
Tongatapu, where wave runups locally exceeded 15 m, are rugged, 
exposed, and mantled by sea cliffs. Only a handful of tourist resorts 
are located there, and these were closed because of the coronavirus 
disease 2019 (COVID-19) pandemic. Similarly, the west coast of 
‘Eua is sparsely populated, limiting human exposure to the 10- to 
12-m waves that beached there. The coast of Tofua is entirely unin- 
habited and rimmed by the highest sea cliffs in the archipelago. The 
main settlements in all of Tonga safely face Nuku'alofa Lagoon, the 
best possible scenario to duck the full brunt of a megatsunami. The 
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many other inhabited islands in the archipelago sit atop the Ha'apai 
platform, also providing ample shallow water to attenuate the 
tsunami to beach with heights <5 m. Again, geography was kind. 
No inhabited islands situate on the western flank of the Ha'apai 
platform where the tsunami were at their fiercest. The normally 
busy resort islands of Ha'ano, Foa, Lifuka, Uolva, and Uiha sit on 
the sheltered eastern flank of the Ha'apai platform. In a lucky twist, 
these resorts were shuttered because of COVID-19. 

Disaster management policies and preparation also contributed 
to the low death toll. After the 2004 Indian Ocean tsunami, coastal 
nations worldwide became more risk alert. The Government of 
Tonga held their first kingdom-wide tsunami evacuation drill in 
November 2017 and performed a comprehensive multidisaster 
risk assessment survey a mere 7 months before the 2022 event. 
These activities, as well as the Polynesian islander's connection 
with their geoheritage (62), are credited as a major component of 
the unexpectedly low number of casualties. The outcomes from 
the HTHH event provide an ideal example of how such efforts are 
critical to mitigate tsunami disaster impacts and to reducing the 
death toll. 

Most volcanically generated megatsunami can only be studied 
through geological or archeological evidence, wherein much of 
the evidence is erased or altered over time. Two well-known exam- 
ples are the 1883 eruption of Krakatoa, wherein 36,000 lives were 
lost (63) and the ~1600 BCE Theran eruption in Santorini, 
wherein the death toll is highly speculative given that only one 
human victim has been reported within well-constrained tsunami 
deposits (64). Both of these events share similar multiple eruption 
sequences similar to observations from HTHH (e.g., four modeled 
for Thera and four described for Krakatoa). HTHH provides the op- 
portunity to observe at scale and in real time through data collected 
at tide and barometric gauges and in the immediate aftermath, pro- 
viding the parameters to better understand those previous events, as 
well as contribute to future disaster preparations. 


MATERIALS AND METHODS 

Topo-bathymetry models 

Key to our tsunami simulation is an accurate topo-bathymetric 
surface. For water depths >25 m, we assembled Global Multi-Reso- 
lution Topography, a global multibeam compilation (65) with a 
spatial resolution of 100 m. Shallower than 25 m, we used the spec- 
trally derived bathymetry of Purkis et al. (20) that uses the depth- 
derivation algorithm of Kerr and Purkis (66) to deliver a spatial res- 
olution of 5 m and is calibrated by abundant sonar soundings from 
the field. 

Coarse-resolution terrestrial DEMs for all islands in the model 
domain were provided by the NASADEM (67) built from the 
ASTER Global Digital Elevation Map (ASTER-GDEM v.3) 
derived from stereo-pairs of satellite scenes from the ASTER 
archive at 30-m spatial resolution (68) combined with the Shuttle 
Radar Topography Mission (SRTM) measurements. Building 
forward from the work of Garvin et al. (50), we mapped those 
islands around which wave runup heights were harvested, using a 
combination of 1- to 3-m compact polarimetric Radarsat Constel- 
lation Mission (RCM) Spotlight C-band SAR (synthetic aperture 
radar) and stereo-pairs of 0.5-m panchromatic Maxar WorldView 
data (WV-02). These terrain models were validated against GEDI 
(69) and ICESat-2 satellite LIDAR. 
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Tsunami simulation 
In the absence of added frictions, the Tsunami Squares method ac- 


celerates fixed squares of fluid at each time step by 
Dvir, t) V P(r, t) 


Dt 


= —gVn(r, t) (2) 


vis fluid velocity (m/s). —g V n(r, t) is the acceleration due to gravity 
g(m/ s*), working on the slope of the water surface Vn(r, t), and — V 
P(r, t)/p,, is the acceleration due to the gradient of explosive over- 
pressure versus time, P(r, t), where P is applied overpressure pres- 
sure (Pa) and ,, is fluid density (kg/mĉ). Data from surface, open- 
air atomic bomb tests (25) fix the decay of peak overpressure with 
distance r as 


1 


W WJ? 
P(r) = 98.368 B + 42.644 E 
r r 


(3) 


Peak overpressure, P, is pounds per square inch (psi = 6895 Pa). r 
is the distance from the blast in kilometers, and W is the explosive 
yield in megatons (1000 kg) of trinitrotoluene (TNT). P(r, t) then 
becomes a function of distance only 


P(r, t) = P(r, t) = P(r)F(W,r, t) (4) 


where F(W, r, t) is a unitless normalized time history of pressure 
that we take as a simple exponential 


E(W,r, t) Se Oe ty (r) (5) 


In Eq. 5, to(r) is the arrival time of the blast at distance r; that is, 
distance in kilometers divided by the speed of sound, 0.34 km/s. The 
decay constant A depends on yield and peak overpressure by 


A(W,r) = 4 max [W"3, 3.87 W!/3/P!/?(p)] (6) 


Units in Eq. 6 are A (in seconds), W (in megatons) and P (in psi). 
Equation 5 is only slightly modified from the work of Newmark 
(26). Because peak overpressure (Eq. 3) goes to infinity as r goes 
to zero, — V P(r, t)/pw in Eq. 2 has to be capped somehow. We 
choose to limit this term to the acceleration due to gravity g. This 
restriction affects only the first few seconds of the blast. 

Last, to account for frictional losses in water <30 m deep, an ad- 
ditional acceleration is added to Eq. 2 


Dy(r, t) V P(r, t) 
Dt 


ca |v(r, t)|v(r, t) 


= —gVn(r, t) An) 


(7) 

w 
where H(r) is water depth and cg is a unitless friction coefficient. 
The simulations used squares of 200-m size and were time 
stepped by 1 s. The behavior of breaking waves and bore formation 
are implicitly captured in the nonlinear formulation of 
Tsunami Squares. 


Validation of the tsunami simulation 

We audited our tsunami simulation across space and through time. 
For the former, we compared simulated versus observed wave 
runups at 118 coastal sites. To assess performance through time, 
we compared water levels recorded by the Queen Salote and Vuna 
Wharf tide gauges on the Nuku'alofa waterfront (Tongatapu) with 
simulated water levels for the duration of the HTHH explosive 
episodes. 
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Of the 118 validation sites, 91 were visited in the 3 months fol- 
lowing the 15 January 2022 HTHH eruption, allowing direct mea- 
surement of wave runups above the tidal stage at the time of the 
tsunami (44). On the basis of sedimentary deposits, erratic clasts, 
floated pumice, and damaged vegetation and infrastructure, two 
wave runup heights (SR, and SR,) were determined for each site 
(fig. S1). These measurements do not represent minimum and 
maximum runup values but rather give a range of potential variabil- 
ity. We also extracted two runup values from the tsunami simula- 
tion; WR, represents the maximum height above sea level where 
flow depth goes to zero (akin to field-measured SR), and WR; is 
like the field-measured SR). 

At a further 27 sites, we assembled pre- and post-tsunami satel- 
lite imagery and estimated runups based on the topography to 
which vegetation had been felled after 15 January 2022. Where 
cloud-free satellite imagery existed before and after the tsunami, 
Maxar WV-02 and WV-03 satellite data were used (Fig. 10). 
These satellites produce images in eight multispectral bands with 
pixel dimensions of 0.5 m by 0.5 m once sharpened using the broad- 
band panchromatic channel. In the rare cases where either the “pre” 
or “post” WorldView optical images were too cloudy, Canadian 
Space Agency's (CSA) RCM-SAR scenes served as replacement 
(Fig. 5). These have pixel dimensions of 1 m by 3 m in the spot- 
light-imaging beam mode (FSL30 at ~45° incidence). Such SAR 
imaging data are uncompromised by clouds and microwave back- 
scatter (70) and easily discern felled vegetation, as well as distur- 
bances related to slopes or changes of meter-scale surface texture. 
In addition, satellite LIDAR measurements of tree heights on the 
southern flanks of Tofua from ISS/GEDI (69) illustrated canopy 
heights in potentially disturbed areas ranging from 5-m to over 
45-m elevation above sea level, with a typical value of about 20 m, 
consistent with highly approximated tree size measurements from 
submeter scale WV images. The 1.5-km section of the south-south- 
west coastline of Tofua where the wave runups were simulated to be 


Elevation (m) £8 


Pre-tsunami 
coastline — 


Post-tsunami 
coastline — 


at their highest was additionally surveyed using a DJI drone. This 
survey delivered a DEM with 8-cm resolution and oblique images 
showing disturbances to the flanks of the island varying from 30 to 
70 m elevation (Fig. 6). Comparison to pansharpened WV-02/03 
images acquired before 15 January 2022 indicate these disturbances 
to postdate the tsunami. For all wave runups derived from remote 
sensing data, we adapted the workflows of Purkis et al. (71) and Wu 
et al. (72) and digitized the most seaward position of standing veg- 
etation in the pre-tsunami image for each site and then repeated this 
exercise for the post-tsunami imagery. We fixed wave runup values 
equal to the highest elevation (as audited from the RCM-SAR 
images and Maxar WorldView DEMs) of vegetation standing 
before 15 January but knocked down in the post-tsunami 
imagery. To estimate variability, we sampled 5-m intervals along 
200-m coastline transects and noted similar heights. The 10th and 
90th percentile of these values form SAT, and SAT,, respectively. 

While a valuable supplement to the field data, estimating wave 
runups from satellite comes with two caveats. First, the spatial res- 
olution of the satellite data limits the fidelity to which coastal change 
can be detected, often because of challenging viewing geometry or 
illumination constraints. Damage to coastal vegetation or deposi- 
tion of erratic clasts at a scale finer than the pixel size of the satellite 
imaging will be missed or ambiguously depicted. Second, minor 
damage to vegetation that might be obvious in the field, is unlikely 
to be detected from orbit, except in special cases via LIDARSs such as 
GEDI, or ICESat-2 (which directly measure vegetation structure). In 
this respect, we consider the 27 satellite-derived runup heights to be 
conservative compared to the 91 field-measured heights. 

Last, to assess the accuracy of the simulation through time, we 
compared water levels recorded at the Queen Salote and Vuna 
Wharfs tide gauges with the levels predicted by the tsunami simu- 
lation. Water levels recorded by the two gauges were detrended to 
remove the tidal signal such that the water level was at 0-m elevation 
at 4:15 UTC. The comparison between the two tide gauges and the 


14 January 2022 Site81 17 January 2022 


14 January 2022 Site82__17 January 2022 


Fig. 10. Pre- and post-tsunami vegetation damage on Tongatapu. (A) Location of three examples from southeastern Tongatapu (site IDs 81, 82, and 83) where 
downed coastal vegetation allows wave runup heights to be estimated from a DEM built from stereo-pairs of WorldView satellite data. (B), (C), and (D) compare the 
position of the vegetated coastline in pre- (14 January 2022) and post-tsunami (17 January) WorldView imagery for these sites and spot measurements of wave 


runup from the DEM. 
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simulation spans the interval 04:15 UTC to 6:00 UTC; the Queen 
Salote Wharf gauge failed at 5:24 UTC. 


Supplementary Materials 
This PDF file includes: 

Fig. S1 

Table $1 

Legends for movie $1 and S2 


Other Supplementary Material for this 
manuscript includes the following: 
Movies S1 and S2 
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